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ABSTRACT 

Context. Internal flows inside gravitationally stable astrophysical objects, such as the Sun, stars 
and compact stars are compressed and extremely subsonic. Such low Mach number flows are 
usually encountered when studying for example dynamo action in stars, planets, the hydro- 
thermodynamics of X-ray bursts on neutron stars and dwarf novae. Treating such flows is nu- 
merically complicated and challenging task 

Aims. We aim to present a robust numerical tool that enables modeling the time-evolution or 
quasi-stationary of stratified low Mach number flows under astrophysical conditions. 
Methods. It is argued that astrophysical low Mach number flows cannot be considered as an 
asymptotic limit of incompressible flows, but rather as highly compressed flows with extremely 
stiff pressure terms. Unlike the pseudo-pressure in incompressible fluids, a Possion-like treat- 
ment for the pressure would smooth unnecessarily the physically induced acoustic perturbations, 
thereby violating the conservation character of the compressible equations. 
Moreover, classical dimensional splitting techniques, such as ADI or Line-Gauss-Seidel methods 
are found to be unsuited for modeling compressible flows with low Mach numbers. 
Results. In this paper we present a nonlinear Newton-type solver that is based on the defect- 
correction iteration procedure and in which the Approximate Factorization Method (AFM) is 
used as a preconditioner This solver is found to be sufficiently robust and is capable of capturing 
stationary solutions for viscous rotating flows with Mach number as small as M ~ 10"'', i.e., 
near the incompressibility limit. 
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1. Introduction 
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Among different energy contents, the gravitational and thermal energies in bound astrophysical 
systems are dominant. The virial theorem states that in the absence of external pressure and surface 
tension the total energy of gravitationally bound system is negative, i.e., 

-f^i-^+2[e,h+£kin]+ySi-^ <0, (1) 
K K 

where ai,/3i are constants less than one and where 

Gravitational energy : £grav — 

Thermal enersy : £th = | L P dvol 

, (2) 
Kinetic energy : £kin — 5 Jy^-'l^fl dvol 

Magnetic energy : £inag - 

where G, M, R, P, Vf, denote respectively the gravitational constant, the magnetic flux, the 
mass and radius of the object, pressure and fluid-velocity, and dvol is an infinitesimal volume- 
element. 

The final stage in the evolution of such gravitationally stable systems is characterized by the fol- 
lowing energy measure: 

Ifigmvl > |£thl » |£kinl, |£mag|. (3) 

In terms of velocities per mass this relation is equivalent to: 

V2>V^»V^,V^, (4) 

where the velocities correspond to the self-gravitating energy (Vg = ^^)^ thermal, fluid and mag- 
netic (Alfven) velocities. 

Therefore, fluid motions in gravitationally stable astrophysical systems are naturally sub-sonic, 
hence the Mach number is relatively low. 

For example, helioseismology measurements have revealed that the Sun oscillates on various fre- 
quencies. In particular, it has been found that the origin of the 5 -minut e oscillations is a self-excited 



1974). This corre- 



sound wave travelling back-and forthwards through the Sun interior (IMusman , 
sponds roughly to the sound speed: 

Vs — ~ 2.3 X 10** cm s"'. (5) 

5 minutes 



Roth et al. 



( 120021) have suggested that internal flows can have a maximum sectorial amplitude of 
about lO^* cm s They argue that a higher velocity would lead to a noticable distortion of the 
rotation rate in the convection zone, hence contradicts observations. In this case, the Mach number 
reads: 

M^^^10-\ (6) 

Consequently, the fluid motions in the Sun is compressible with extremely low Mach numbers. 
Similarly, in the case of neutron stars, the te mperature of the s uperfluid ranges between 10^ up to 



5 X 10 K, depending on the crust heat source jVan Riperi,ll991h . The superfluid velocity relative to 
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2003h . 



coordinates rotating with angular velocity Qns can reach Vhd ~ 10'* - 10'' cm s ' (I Jones , 
Thus, the ratio of the sound crossing time to the hydrodynamical time scale reads: 

I^^(Y^]\m'^10-\ (7) 

where V^[- dP/dp - (p - \p)l{p + p)] corresponds to the the sound speed squared, which is 
roughly 10% the speed of light, depending on the equation of state. 

The flows in these two extreme astrophysical objects indicate that numerical solvers should be ro- 
bust enough to deal with extremely low Mach number flows. Such flow-conditions are encountered 
when trying to model the origin of the solar dynamo or the thermonuclear ignition of hydrogen 
rich matter on the s urface of neutron stars, considered to be responsi ble for Type-I X-ra y bursts 



dFisker et al. 



20051) or for novae eruption in the case of white dwarfs dCamenzind 



2007). 



2. Compressible versus weakly and strongly incompressible flows 

While the equations describing compressible and incompressible flows are apparently similar, the 
underlying physics and the corresponding numerical treatments are fundamentally difi'erent. 
In general, compressible flows are made of plasmas. The internal macroscopic motions may be- 
come either supersonic or extremely subsonic. Incompressible flows however, are generally made 
of liquid, so that a further compression would not lead to a noticeable change of their density. The 
transition from gas phase into fluid phase mostly does not occur via smooth change of the equation 
of state. For example, a high pressure acting onto a container of hot water vapor cannot be asymp- 
totically extended to describe the pressure in normal water fluid. Therefore, from the astrophysical 
point of view, weakly incompressible flows can be viewed as strongly compressed plasmas, in 
which the macroscopic velocities are relatively small compared to the sound velocity. 
To clarify these differences, we write the set of hydrodynamical equations in non-dimensional form 
using the scaling variables listed in Table ([1]). 

Scaling variables neutron star(interior) 



L Length ~ 10'^ cm 

p Density ~ lO''' g cm"^ 

T Temperature ~ 10^ K 

T Pressure ~ 10^' dyn cm"^ 



V Velocity ~ 10'' cms"' 

B Magnetic Fields ~ 10' G 
M Mass ~ Mq 

Table 1. Scaling variables for non-dimensioning the hydrodynamical equations. 



The set of hydrodynamical equations describing compressible plasmas in conservative form reads: 

- Continuity equation: 

^ + v-py = o, (8) 

ot 

- The momentum equations: 

dpV 1 p /Mmag\^ 1 

dt ^ Af Fi^ \ M Re 
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where cr{= rjiW + (W)^) - |?7(V • V)!), J] = pv and , are the Reynolds stress tensor, the 
dynamical viscosity coefficient and the gradient of the potential energy, respectively. 

- The total energy equation: 

^ + V • (6 + p)y = pVY • y + i^f V • (Vcr) + ^ V . (vtVD, (10) 

at Fr Re Pe 

where fi = p(e + ^V^) and vt is the heat diffiision coefficient. 

We may simplify the total energy equation by separating the internal energy from the mechani- 
cal energy and assuming a perfect conservation of the latter. Hence, we are left with an equation 
that describes the time-evolution of the internal energy: 

_ + v.£'^y = -(7-i)£''v.y + (r-i){l— J t + — v • (vtYD}, (ii) 

where T(= 77IV • V\^) is the dissipation function. 

- Magnetic equation 

The magnetic induction equation, taking into account transport and diffusion in non- 
dimensional form reads: 
dB 1 A1™8 

^=V>«^x«-r;^(^)Vx«>- (12) 



Name Symbol Definition 



Reynolds number 


Re 


VL/v 


Mach number 


M 


V/Vs 


Reynolds number (magnetic) 


J^ginag 


V L/v™"* 


Mach number (magnetic) 




Va/Vs 


Prantl number 


Pr 


v/vt 


Froude number 


Fr 




Peclet number 


Pe 


RePr 



Table 2. Nondimensional numbers. In this list, the additional parameters v, V^^, vj and Vg cor- 
respond to hydrodynamical viscosity coefficient, magnetic diffusivity coefficient, heat diffusion 
coefficient and to the effective velocity of the potential energy *P(i.e., = VP, respectively. 



On the other hand, incompressible flows are described through the following set of equations: 
V-y = 0. (13) 

M^p [V, + (y • V)y] = -VP + — rp + -—V • a (14) 

Fr^ Re 

f .v.ry = (r-i){(^)'T.lv.(v,vr)}, (i5) 

Despite the apparent similarity, the pressure in compressible flows has different physical meaning; 
in the compressible case the equation of energy influences the momentum equation through the 
equation of state, whereas in the incompressible case, the pressure is just a lagrangian multiplier 
with no direct physical meaning. The set of equations of incompressible flows is characterized by 
the following two features: 

1. The velocity field must not only evolve as described by the momentum equations, it should 
fulfill the divergence-free condition also. 
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2. there is no direct equation that describes the time-evolution of the pressure. 

Therefore, we may use the pressure in the momentum equations to form an equation that enforces 
the flow to move in such a manner that the divergence-free condition is always fulfilled, indepen- 
dent of the constitutive nature of the flow. 

This can be achieved by taking the divergence of the momentum equation above: 



M'pW, + iV ■ V)y] - -^p V>P - -^V ■ cr 
which is equivalent to the following compact form: 
AP ^RHS. 



= -V ■ Vf = -AP, 



(16) 



(17) 



The right hand side (RHS) contains the divergence of the other terms of the momentum equation. 



The strategy of turning VP in the momentum equations into a Possion equation to achieve 
divergence-free motions is the basis of different variants of the so called projection method, e.g., 
the "Semi-Implicit Method for Pressure-Link ed Equation" (SIMPLE) and "Pressure-Implicit with 
Splitting Operator" (PISO, see iBartonll 19981 for further details). Similarly, the projection method 



can be applied also to the induction equation in magnetohydrodynamics. Here the induction equa- 
tion is modified to include the gradient of a scalar function as follows: 



— = Vx<yxB 
dt 



■> + V0. 

Taking the divergence of this equation, we obtain: 
A© = V ■ B. 



(18) 



(19) 



Therefore, this method violates the conservation of the magnetic flux. The V0 is actually a source 
function for generating or annihilating magnetic flux, so that generating of magnetic monopoles 
from zero magnetic flux cannot be excluded. Specifically, if V0 = const., then A0 = 0, which 
means that a constant pumping of magnet flux due to numerical errors cannot be eliminated by 
applying a Poisson like-operator 

The matrix form of the projection method applied to the momentum Equation ( fT4l i and to the 
Possion Equation (T3[ is as follows: 

yn+l 
pn+1 



J G 

G* 



RHS 




(20) 



where the coefficient matrices J - dL^n/dV, G - dL^ldP, G* = dL^/dV, and "n+V corresponds 
to the values at the new time level. Applying LU-decomposition, the matrix Equation ( |20] | can be 
re-written as: 



J 

G* -GT^G 



I J-^G 
/ 



yn+l 
pn+1 



RHS 




(21) 



where I denotes the identity matrix. This equation is solved in two steps: 



J 

G* 





-G'J-^G 



V* 
P* 



RHS 




II. 



I 



yn+l 
nn+l 



V* 
P* 



(22) 
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In general the inversion of the Jacobian, J, is difiicuh and costly, and it is therefore suggested to 
replace it by the preconditioning A. In this case the above-mentioned two step solution procedure 
should be reformulated and solved using the defect-correction iteration procedure: 



II. 



A 

G* -G'A-^G 
I A-^G 
/ 



6V* 
6P* 



where (5P* = P* - P", (5P" 



6V* 
6P* 



(23) 



P* andd = RHS = JV*. 



Based on an extension of the projection method, the Multiple Press ure Variable m ethod for mod- 
eling weakly incompressible flows has been suggested (MPV, see 



Munz . 



20031) . Following this 

scenario, it is argued that in the low Mach number regime (M « 1), the variables can be ex- 
panded in the following manner: 



(24) 



However, contrary to the flow conditions in the Sun or violent fluid motions on the surface of com- 
pact objects powered by thermonuclear flashes, the MPV method requires to set the first leading 
terms in the expansion of the pressure to zero in order to assure matching of the solutions in the 
asymptotic limit. Moreover, the MPV expansion requires that the Mach number be sufficiently 
small in order to get an adequate asymptotic limit. 

Almegren et al. 2006 studied the time-evolution of an injected heat bubble in the atmosphere of a 
neutron star using different types of numerical approximations aimed at properly treating incom- 
pressible flows in the weak regime. They find that their suggested low Mach number approximation 
performs relatively well compared to pure incompressible or anelastic approximations. The strat- 
egy relies on finding an appropriate function /3 of the initial conditions such that V ■ /3V = is 
always fulfilled. The leading terms of the pressure is set to oppose gravity, so that the core remains 
in hydrostatic equilibrium. 

We note however that the violent flow-motions associated with Type-I X-ray bursts may not nec- 
essary be of low Mach number type, so that finding the appropriate /3 may turn into a difficult 
analytical task. 



3. Highly compressed low Mach flows: an iterative non-linear preconditioned 
Newton solver 



Assume we are given a two-dimensional nonlinear vector equation of the form: 
dq dF(q) dG(q) 



dt dx 



dy 



(25) 



where q, F, G, f denote the vector of variables, their momentum flux both in x and y-directions and 
a source function, respectively. 

Define the residual d(q) and look for the vector q, such that d(q) - 0, i.e.. 



d(q) = / ■ 



dq dF(q) dG(q) 
dt dx dy 



0. 



(26) 
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.... 

.ii"ip' "i.",', "li,, 
■IP" 



real Jacobian 
= J 



p''"iS:!l'v:mji 

. ji _.■ ^ . X 



Approximate Jacobian 
= J 



■ H % 



ID: low resolution 



Fig. 1. The real Jacobian matrix (J, left) and an approximate Jacobian matrix (J, right). In J we have displayed 
the entries (blocks) that correspond to low and high spacial resolution in one and two dimensions. 



In the case of a single one-dimensional nonlinear function Tix) — 0, the zeros can be found using 
the Newton iteration method: 



(27) 



where ;F(x') = ^lx=x' "i" denotes the iteration number. When applying this approach to a 
general system of equations such as Eq. dZST l. we have then to perform the following replacements: 



X ^ q, 
r(x) d{q) 
T ^ J, 



r'd 



(28) 



where J is the Jacobian matrix defined as: 7 = |^ 



dq- 



Defining - q^^^ 
J;u = d, 



■ q^, we may re-write Eq. 



as: 



(29) 



(30) 



where "d" is calculated using arbitrary high spatial and temporal accuracies. 
The matrix Equation ( |29] l is said to be: 

I Linear: if d = d(q") 
I Otherwise : if d = d(q'). 

While in the first case, one need to invert the Jacobian once per time step, in the second case how- 
ever, several iterations per time step might be required to recover the nonlinearity of the solution. 
The calculation may become prohibitively expensive, if the Jacobian to be inverted corresponds to 
a system of equation in multi-dimensions to be solved with high spatial and temporal accuracies. 
The idea of preconditioning is to calculate the defect "d" as proposed by the physical problem (e.g., 
with very high resolution), whereas the Jacobian is replaced then by an approximate matrix A of 
the following properties: 

- A is easier to invert than J, 

- A and J are similar and share the same spectral properties. 

While the first property is easy to fulfill, the second one is in general an effort-demanding issue. It 
states that the preconditioning A should diff'er only slightly from the Jacobian if trivial replacement 
to be avoided, therefore, given the matrix A, the solution procedure would run as follows: 
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1. Compute the defect "d". 

2. Use the matrix equation: Aju - d.to solve foryU. 

3. Update: q'^^ - + jj, re-calculate "d" and A, respectively. 

4. The procedure (2) and (3) should be repeated until max(|(i|) is smaller than a number e, where 
the maximum runs over all the elements of "d". 

The fundamental question to be addressed still is: how to construct a robust preconditioner A that 
is capable of modeling low-Mach number flows efficiently, but still easy to invert? 
In this construction, two essential constrains should be taken into account: 

- A conservative first order spatial discretization of the Navier-Stokes equations generally yields 
a Jacobian matrix of penta-diagonal block form as depicted in Figure ([TJ. 

- The gradients of the thermal pressure are dominant, so that the pressure-related terms must be 
treated simultaneously. 

In order to clarify these points, we re-write the matrix Equation ( |29] l at an arbitrary grid point (j,k) 
in the following block form: 

+5:>M,k +£>j7Vj,k +5-kyUj+i,k = ^/j,k (31) 

where S_'^'^, S denote the sub and super-diagonal block matrices and - II5t + + the 
diagonal block matrices, respectively. 

While this block structure is best suited for using the one-colored or multi-colored line Gauss- 
Seidel iterative method, test calculations have shown however, that these iterative methods may 
stagnate or they may even diverge if the flow is of low mach number type. The reason for this 
behaviour is that most iterative methods rely either on partial updating of the variables or on the 
dimensional splitting. These, however, are considered to be inefficient methods or they may even 
stagnate, if the system of equations to be solved are of elliptic type, such as the Possion equations. 
One way to take the multi-dimensional variations of the pressure all at one time is to spatial- 
factorize the Jacobian into sub-matrices, such that the resulting multiplication results in a good 
approximation of the original Jacobian, i.e., 

J M n^AiA^. (32) 

The advantages of this method is that the gradients of the pressure in all direction are incorporated 
in the matrix A. Furthermore, the right hand side, i.e., the defect is updated only after all matrices 
Am are fully inverted. We note that this preconditioner can be used even as a direct solver as long 
as time-dependent solutions are concerned. To clarify this procedure, we rewrite Equation ( |25T l in 
the finite space as follows: 

So A^F"+i A„G"+i 

— + + = 0. (33) 

5t Ax Ay 

where (5q - q"^' - q"(//) and Ax,y are space difference operators. 

We may expand F°^' and G"^' around their values at the old time levels as follows: 
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Equivalently, 

p„+l ^ pn ^ ^«^^ ^ ^(^^^2^ 

G"+' = G" + B"Sq + OiStf, 

Substituting these expressions into Equation (|33] l. we obtain: 



- + LxA" + LyB" + 0(5tY 



6q = UF" + LyG", 



(34) 



(35) 



where Lx, Ly denote the differential operators in x and y-directions, respectively. We may replace 
Equation ( [35] ) by the following approximation: 

" I 



6i 



+ LxA" 



[l + 6t LyB"] 5q = LxF" + LyG", 



(36) 



This replacement induces an error which is proportional to : dtL^Ly + 0{6t)^. This error may 
diminish for steady conserved fluxes, but may diverge for time-dependent solutions if the time 
steps are large. The latter disadvantage is relaxed by the physical consistency requirement, that 
small time steps are to be used if the the sought solutions are time-dependent. 
The matrix Equation ( l36b can be re-written in the following compact form: 



Ax Ayju = d. 



(37) 



where Ax = + UA" and Ay = + LyB". 



Applying this factorization within a non-linear iterative solution procedure, the solution procedure 
would run as follows: 

1 . Compute the defect d - d(q\ q") at each grid point using the best available spatial and temporal 
accuracies. 

2. Solve the matrix equation: Ay^dq* = d, to obtain 6q*. 

3. Solve the matrix equation: Aydq - 6q* to obtain 6q. 

4. Update q: - q^ + dq and subsequently the defect d. 

5. Perform a convergence check to verify if max(|c/|) < e. If not, then repeat the algorithmic steps 
1-4 repeatedly. 




Fig. 2. Two concentric spheres: the inner sphere has a radius Ri and rotates with angular velocity Oi whereas 
the outer one has the radius Ri and rotates with 
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4. Taylor-Couette flows between two concentric spheres 

Large scale motions of gas in stellar spherical shells are controlled by the imbalance of energies, 
namely between the potential, thermal, rotational and magnetic energies. It is generally accepted 
that rotation deforms surfaces of constant pressure, but has only indirect influence on surfaces of 
constant te mperatu r es. Th e resulting barocUnicity is unbalanced and derives large scale meridional 
circulation (ISweeti Il950l) . On the local scale, these flows are in general convectively unstable. 



hence governed by convective turbulence. Such combined motions are observationally evident in 
the solar convective zone. 

In the laboratory, spherical Couette flows between two rotating spheres are considered to be similar 
to rotating stellar envelopes. In the case of fast rotation, the flow is a combination of primary az 



imuthal rotations and a secondary meridional circulation induced by Ekman pumping ([Greenspan, 



19681) . Here the flow is controlled by two parameters: The Reynolds number and the gap width 



between the two spheres. The Reynolds number for this configuration is defined as: 

|An|Ri|AR| 1^2- nil Ri \R2-Ri\ 
Re , = , (38) 

V V 

where R12, fii,2 and v are the inner and outer radii, the angular frequency of the inner and outer 
sphere and the viscosity coefficient, respectively (Figure|2]i. 

The number of rotationally-induced fluid vortices and transition to turbulence in Couette flows 
depend on how large the Reynolds number is as well on the width S of the gap between the two 
spheres. For example, for Re > 460 and 6 - 1 .006 Couette flows have been verified to become 



turbulent (IGertsenshtein et al. 



2OOII) . 



In applying our solver to Couette flows between two-concentric spheres, the following parameters 
have been used: Qi = 3, Q2 = 0, Ri - I, 7?i = 1.25 and a viscosity coefficient v = 0.005 (see 
Figure |3j. 

The domain of calculation is limited to the first quadrant [1 < /? < 1.25] x[Q < 6 < n/l]. Along the 
equator and polar axis reflecting boundary conditions have been imposed, whereas the outer and 
inner boundaries are set to be rigid with zero material flux across them. 

In order to test the capability of the solver to deal with low Mach number flows, we have system- 
atically increased the temperature from 10 up to one million, which corresponds to a reduction of 
the Mach number by at least two orders of magnitude. 

Our results, partially displayed in Figure (|3]l, indicate that AFM as a preconditiong in combination 
with the defect-correction Newton iteration is indeed capable of modelling weakly incompressible 
flows down to Mach number M ~ 10^^. 



5. Summary 

In this paper we have shown that weakly incompressible flows in astrophysics are actually highly 
compressible low Mach number flows in which the pressure plays a vital role in dictating the fluid 
motions. Therefore, these flows can be well-treated using a robust compressible numerical solver 
in combination with a sophisticated treatment of the pressure terms. 

However, as the sound wave crossing time in low Mach number flows is extremely short relative to 

the hydrodynamical time scale, we have concluded that time-explicit methods are not suited. 

On the other hand, projection methods based on Possion-like solvers for the pressure violates the 
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Fig. 3. Taylor-Couette flows: the inner sphere has a radius Ri = 1 and rotates with £li = 3, whereas the outer 
sphere has a radius Ri = 1.2 and fli = 0. The flow has the constant viscosity coefficient a = 0.005. The 
Mach number is set to decrease systematically by increasing the temperature from T = 10^ up to T = 10*. 

This corresponds to a reduction of the Mach number by two orders of magnitude. The right panel shows the 
corresponding time-evolution of the time step size in units of Courant number. 

conservative formulation of the hydrodynamical equations, namely the entropy principle and there- 
fore the monotonicity character of the scheme. 

Our conclusion is that the set of hydrodynamical equations describing the time-evolution of com- 
pressible low Mach number flows should be solved using an impUcit robust solver. Our numerical 
studies show that nonhnear Newton-type solvers in combination with the defect-correction itera- 
tion procedure and using the Approximate Factorization Method (AFM) as preconditioning is best 
suited for treating such flows. Unlike the classical non-direct methods that rely on dimensionjil 
splitting and/or partial updating, e.g., ADI and Line-Gauss-Seidel, the AFM is based on factoriz- 
ing the Jacobian matrix and subsequently updating the variable in all directions simultaneously. 
Similar to the treatments of the pseudo-pressure in incompressible flows, the pressure gradients 
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in low Mach number flows do not accept dimensional splitting even when it is applied within the 
preconditioning only. 
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